mydir.scripts <- "./scripts"
mydir.output <- "./output"
mydir.data <- "./data"
mydir.output.perGene <- mydir.output
source(paste0(mydir.scripts, "/functions_wilcoxon_t_logistic_fdr.R"))
The LLS data was produced as part of the BBMRI project. Data can be requested for registered users via their website.
The GEUVADIS data has been pre-processed again in the same way as the LLS data, and analyses were rerun. The original GEUVADIS data can be obtained for example from the EBI.
The script used to install the spliceQTL package, obtain SNP and exon-level data, and produce objects needed, can be found on Armin Rauschenberger’s own [github](https://github.com/rauschenberger/spliceQTL/blob/master/vignettes/code.Rmd] page. Lines 1-131 refer to installing, preprocessing and running analyses. In the remainder results are explored. Some of the results in the article are included below.
Researchers willing to perform these analyses for other datasets are encouraged to use the chunks in that report.
Below we reproduce some of the graphs presented in our article, loading results obtained.
The spliceQTL test was run in the same way as done for the analysis of the GEUVADIS data (with the original pre-processing). Here we load the results of that analysis.
vchrs <- 1:22
mtable <- matrix(0, nrow = length(vchrs), ncol = 4)
colnames(mtable) <- c("none", "onlyGEU", "onlyLLS", "both")
rownames(mtable) <- vchrs
mtable <- data.frame(mtable)
table.g <- table.l <- table.go <- table.lo <- NULL
mlist <- list()
for(chr in vchrs)
{
load(paste0(mydir.data, "/pvalNew.Geuvadis.chr", chr, ".RData")) # pvalue
a <- pvalue; pvalue <- NA
a <- data.frame(a, chr = rep(chr, length(a)))
table.g <- rbind(table.g, a)
load(paste0(mydir.data, "/pvalNew.LLS.chr", chr, ".RData")) # pvalue
b <- pvalue; pvalue <- NA
b <- data.frame(b, chr = rep(chr, length(b)))
table.l <- rbind(table.l, b)
# filter tests (trial)
names <- intersect(rownames(a),rownames(b))
a <- a[names, ]
b <- b[names, ]
# If desired, compute correlations between results of the 2 datasets
if(FALSE){
mtable$cor[ chr ] <- round(stats::cor(-log(a[, 1]), -log(b[, 1]), method="pearson"), digits=2)
mtable$cor.test[ chr ] <- signif(stats::cor.test(-log(a[, 1]), -log(b[, 1]), method="pearson")$p.value, digits=2)
}
# significance
a <- stats::p.adjust(a[, 1]) < 0.05
b <- stats::p.adjust(b[, 1]) < 0.05
mtable$none[ rownames(mtable) == chr ] <- sum(!a & !b, na.rm = TRUE)
mtable$onlyGEU[ rownames(mtable) == chr ] <- sum(a & !b, na.rm = TRUE)
mtable$onlyLLS[ rownames(mtable) == chr ] <- sum(!a & b, na.rm = TRUE)
mtable$both[ rownames(mtable) == chr ] <- sum(a & b, na.rm = TRUE)
}
save(table.g, file=paste0(mydir.output, "/table_allpvalues_GEUV.RData"))
save(table.l, file=paste0(mydir.output, "/table_allpvalues_LLS.RData"))
# Tables with probes in the overlap are saved in the next chunk,
# after computing the FDR
save(mtable, file=paste0(mydir.output, "/table_cross_GEUV_LLS.RData"))
load(paste0(mydir.output, "/table_allpvalues_GEUV.RData")) # table.g
load(paste0(mydir.output, "/table_allpvalues_LLS.RData")) # table.l
load(paste0(mydir.output, "table_cross_GEUV_LLS.RData")) # mtable
mcut <- 0.05
dname <- "GEUVADIS"
mytable <- table.g
fdr <- NULL
for(chr in vchrs)
{
mname <- paste(dname, "Chr", chr)
vfdr <- adj.bhfdr(mytable[mytable$chr == chr, 1], myvar = mname, cutoff = mcut)
fdr <- c(fdr, vfdr)
}
mytable$fdr <- fdr
table.g <- mytable
###
dname <- "LLS"
mytable <- table.l
fdr <- NULL
for(chr in vchrs)
{
mname <- paste(dname, "Chr", chr)
vfdr <- adj.bhfdr(mytable[mytable$chr == chr, 1], myvar = mname, cutoff = mcut)
fdr <- c(fdr, vfdr)
}
mytable$fdr <- fdr
table.l <- mytable
names.o <- intersect(rownames(table.l),rownames(table.g))
table.lo <- table.l[names.o, ]
table.go <- table.g[names.o, ]
save(table.go, file=paste0(mydir.output, "/table_overlappvalues_GEUV.RData"))
save(table.lo, file=paste0(mydir.output, "/table_overlappvalues_LLS.RData"))
In the code chunk below plots of the FDRs per chromosome and dataset are produced. These make use of all results per dataset, not the results from the overlapping genes. Graphs are saved into a pdf file.
myn <- "GEUV"
myd <- table.g
pdf(paste0(mydir.output, "/plot_fdr_allChrs_", myn, ".pdf"), width = 8, height = 12)
par(mfrow = c(6, 4), mar = c(5, 4, 0.5, 0.5))
for(chr in vchrs)
{
mf <- myd$fdr[myd$chr == chr]
plot(sort(mf), type = "l", col = "blue", xlab = "genes", ylab = "FDR")
mfs <- mf[mf <= 0.05]
lines(sort(mfs), col = "red", lwd = 2)
segments(0, 0.05, length(mf), 0.05, lty = "dashed", col = "gray")
text(0, 0.7, labels = paste(length(mfs), "FDR<0.05"), pos = 4, cex = .8)
text(0, 0.9, labels = paste("chr", chr), pos = 4, cex = .8)
}
dev.off()
myn <- "LLS"
myd <- table.l
pdf(paste0(mydir.output, "/plot_fdr_allChrs_", myn, ".pdf"), width = 8, height = 12)
par(mfrow = c(6, 4), mar = c(5, 4, 0.5, 0.5))
for(chr in vchrs)
{
mf <- myd$fdr[myd$chr == chr]
plot(sort(mf), type = "l", col = "blue", xlab = "genes", ylab = "FDR")
mfs <- mf[mf <= 0.05]
lines(sort(mfs), col = "red", lwd = 2)
segments(0, 0.05, length(mf), 0.05, lty = "dashed", col = "gray")
text(0, 0.7, labels = paste(length(mfs), "FDR<0.05"), pos = 4, cex = .8)
text(0, 0.9, labels = paste("chr", chr), pos = 4, cex = .8)
}
dev.off()
Tables of significant results from the overlapping genes. First without multiple testing correction, and testing for association between the two test results. This is supplementary figure 10.
table.l <- table.l[rownames(table.l) %in% rownames(table.g), ]
table.g <- table.g[rownames(table.g) %in% rownames(table.l), ]
#all(rownames(table.g) == rownames(table.l))
dcols <- densCols(-log(table.g$a), -log(table.l$b))
plot(-log(table.g$a), -log(table.l$b), pch = 20, col = dcols,
xlab = "Geuvadis", ylab = "LLS", main = "SpliceQTL p-values per gene")
abline(h = 3, lty = "dashed", col = "gray", lwd = 2)
abline(v = 3, lty = "dashed", col = "gray", lwd = 2)
Tables, no multiple testing correction. Different thresholds.
vth <- c(0.001, 0.005, 0.01)
for(th in vth)
{
myt <- table(table.g$a <= th, table.l$b <= th)
print(myt)
print(chisq.test(myt)) # test of marginal independence, not quite what we want
print(mcnemar.test(myt)) #test of symmetry
}
##
## FALSE TRUE
## FALSE 14915 1441
## TRUE 421 603
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: myt
## X-squared = 2323.8, df = 1, p-value < 2.2e-16
##
##
## McNemar's Chi-squared test with continuity correction
##
## data: myt
## McNemar's chi-squared = 557.66, df = 1, p-value < 2.2e-16
##
##
## FALSE TRUE
## FALSE 14168 1811
## TRUE 590 811
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: myt
## X-squared = 2175.5, df = 1, p-value < 2.2e-16
##
##
## McNemar's Chi-squared test with continuity correction
##
## data: myt
## McNemar's chi-squared = 619.91, df = 1, p-value < 2.2e-16
##
##
## FALSE TRUE
## FALSE 13213 2249
## TRUE 830 1088
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: myt
## X-squared = 1954.2, df = 1, p-value < 2.2e-16
##
##
## McNemar's Chi-squared test with continuity correction
##
## data: myt
## McNemar's chi-squared = 653.04, df = 1, p-value < 2.2e-16
vth <- c(0.001, 0.005, 0.01)
for(th in vth)
{
myt <- table(table.g$a <= th, table.l$b <= th)
print(myt)
# print(chisq.test(myt)) # test of marginal independence, not quite what we want
# print(mcnemar.test(myt)) #test of symmetry
myt <- 100*myt/sum(myt)
print(chisq.test(myt)) # test of marginal independence, not quite what we want
print(mcnemar.test(myt)) #test of symmetry
}
##
## FALSE TRUE
## FALSE 14915 1441
## TRUE 421 603
## Warning in chisq.test(myt): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: myt
## X-squared = 9.0074, df = 1, p-value = 0.002689
##
##
## McNemar's Chi-squared test with continuity correction
##
## data: myt
## McNemar's chi-squared = 2.2127, df = 1, p-value = 0.1369
##
##
## FALSE TRUE
## FALSE 14168 1811
## TRUE 590 811
## Warning in chisq.test(myt): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: myt
## X-squared = 9.1674, df = 1, p-value = 0.002464
##
##
## McNemar's Chi-squared test with continuity correction
##
## data: myt
## McNemar's chi-squared = 2.628, df = 1, p-value = 0.105
##
##
## FALSE TRUE
## FALSE 13213 2249
## TRUE 830 1088
## Warning in chisq.test(myt): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: myt
## X-squared = 8.7048, df = 1, p-value = 0.003174
##
##
## McNemar's Chi-squared test with continuity correction
##
## data: myt
## McNemar's chi-squared = 2.8975, df = 1, p-value = 0.08872
The correlation between the two test results (on the -log scale) is 0.43081 for Pearson and 0.2708648 for Spearman.
Now with multiple testing correction.
table(table.go[, "fdr"] <= 0.05, table.lo[, "fdr"] <= 0.05)
##
## FALSE TRUE
## FALSE 13954 2137
## TRUE 486 803
for(chr in vchrs)
{
print(chr)
print(table(table.go[table.go$chr == chr, "fdr"] <= 0.05, table.lo[table.lo$chr == chr, "fdr"] <= 0.05))
}
## [1] 1
##
## FALSE TRUE
## FALSE 1526 194
## TRUE 39 76
## [1] 2
##
## FALSE TRUE
## FALSE 980 123
## TRUE 40 61
## [1] 3
##
## FALSE TRUE
## FALSE 773 88
## TRUE 27 38
## [1] 4
##
## FALSE TRUE
## FALSE 523 71
## TRUE 22 32
## [1] 5
##
## FALSE TRUE
## FALSE 661 102
## TRUE 28 34
## [1] 6
##
## FALSE TRUE
## FALSE 741 86
## TRUE 33 37
## [1] 7
##
## FALSE TRUE
## FALSE 669 93
## TRUE 31 42
## [1] 8
##
## FALSE TRUE
## FALSE 518 80
## TRUE 16 28
## [1] 9
##
## FALSE TRUE
## FALSE 564 74
## TRUE 14 17
## [1] 10
##
## FALSE TRUE
## FALSE 594 78
## TRUE 13 26
## [1] 11
##
## FALSE TRUE
## FALSE 879 135
## TRUE 26 44
## [1] 12
##
## FALSE TRUE
## FALSE 820 93
## TRUE 36 49
## [1] 13
##
## FALSE TRUE
## FALSE 265 22
## TRUE 7 12
## [1] 14
##
## FALSE TRUE
## FALSE 486 65
## TRUE 13 28
## [1] 15
##
## FALSE TRUE
## FALSE 403 83
## TRUE 23 33
## [1] 16
##
## FALSE TRUE
## FALSE 603 141
## TRUE 14 36
## [1] 17
##
## FALSE TRUE
## FALSE 801 191
## TRUE 27 53
## [1] 18
##
## FALSE TRUE
## FALSE 222 27
## TRUE 10 14
## [1] 19
##
## FALSE TRUE
## FALSE 1050 226
## TRUE 39 81
## [1] 20
##
## FALSE TRUE
## FALSE 436 62
## TRUE 11 15
## [1] 21
##
## FALSE TRUE
## FALSE 134 31
## TRUE 4 11
## [1] 22
##
## FALSE TRUE
## FALSE 306 72
## TRUE 13 36
Now Manhattan plots are produced, saved into pdfs for 5 consecutive chromosomes at a time. These correspond to supplementary figures 11–15.
# This uses the overlap
dname <- "GEUVADIS and LLS"
pdf(paste0(mydir.output, "/plot_manhattan_both_1.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[1:5])
{
mname <- paste(dname, "Chr", chr)
a <- table.go[table.go$chr == chr, "fdr"]
b <- table.lo[table.lo$chr == chr, "fdr"]
graphics::plot.new()
graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
graphics::box()
graphics::axis(side=1)
graphics::axis(side=2)
graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
graphics::abline(h=0, col="grey")
title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()
pdf(paste0(mydir.output, "/plot_manhattan_both_2.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[6:10])
{
mname <- paste(dname, "Chr", chr)
a <- table.go[table.go$chr == chr, "fdr"]
b <- table.lo[table.lo$chr == chr, "fdr"]
graphics::plot.new()
graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
graphics::box()
graphics::axis(side=1)
graphics::axis(side=2)
graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
graphics::abline(h=0, col="grey")
title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()
pdf(paste0(mydir.output, "/plot_manhattan_both_3.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[11:15])
{
mname <- paste(dname, "Chr", chr)
a <- table.go[table.go$chr == chr, "fdr"]
b <- table.lo[table.lo$chr == chr, "fdr"]
graphics::plot.new()
graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
graphics::box()
graphics::axis(side=1)
graphics::axis(side=2)
graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
graphics::abline(h=0, col="grey")
title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()
pdf(paste0(mydir.output, "/plot_manhattan_both_4.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[16:20])
{
mname <- paste(dname, "Chr", chr)
a <- table.go[table.go$chr == chr, "fdr"]
b <- table.lo[table.lo$chr == chr, "fdr"]
graphics::plot.new()
graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
graphics::box()
graphics::axis(side=1)
graphics::axis(side=2)
graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
graphics::abline(h=0, col="grey")
title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()
pdf(paste0(mydir.output, "/plot_manhattan_both_5.pdf"), width = 8, height = 12)
par(mfrow = c(5, 1), mar = c(2, 3, 3, 1))
for(chr in vchrs[21:22])
{
mname <- paste(dname, "Chr", chr)
a <- table.go[table.go$chr == chr, "fdr"]
b <- table.lo[table.lo$chr == chr, "fdr"]
graphics::plot.new()
graphics::plot.window(xlim=c(0.5, length(a)+0.5), ylim=c(-1,1)*max(-log10(c(a, b))))
graphics::box()
graphics::axis(side=1)
graphics::axis(side=2)
graphics::segments(x0=seq_along(a), y0=0, y1=-log10(a), col="darkred")
graphics::segments(x0=seq_along(b), y0=0, y1=log10(b), col="darkblue")
graphics::abline(h=0, col="grey")
title(paste("Chr", chr, "FDR - Geuvadis (red) and LLS (blue)"))
}
dev.off()
Scatterplots of raw p-values, as well as FDR-corrected p-values, per chromosome are produced below. Results are saved into pdf files.
# pvalues
dname <- "GEUVADIS and LLS"
pdf(paste0(mydir.output, "/plot_pvalues_GEUV_LLS_perChr.pdf"), width = 10, height = 5)
par(mfrow = c(1, 2))
for(chr in vchrs)
{
mname <- paste(dname, "Chr", chr)
a <- -log10(table.go[table.go$chr == chr, 1])
b <- -log10(table.lo[table.lo$chr == chr, 1])
mylims <- range(c(a, b))
dcols <- densCols(a, b)
plot(a, b, col = dcols, pch = 20, main = mname, xlab = "Geuvadis -log10 pvalue", ylab = "LLS -log10 pvalue")
segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)
}
dev.off()
# fdr
pdf(paste0(mydir.output, "/plot_fdrs_GEUV_LLS_perChr.pdf"), width = 10, height = 5)
for(chr in vchrs)
{
mname <- paste(dname, "Chr", chr)
a <- -log10(table.go[table.go$chr == chr, "fdr"])
b <- -log10(table.lo[table.lo$chr == chr, "fdr"])
mylims <- range(c(a, b))
dcols <- densCols(a, b)
plot(a, b, col = dcols, pch = 20, main = mname, xlab = "Geuvadis -log10 FDR", ylab = "LLS -log10 FDR")
segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)
}
dev.off()
Scatterplots of the FDR-corrected p-values of both datasets are also produced, and saved into a pdf file.
# fdr
pdf(paste0(mydir.output, "/plot_fdrs_bothDatasets.pdf"), width=8, height = 8)
a <- -log10(table.go[, "fdr"])
b <- -log10(table.lo[, "fdr"])
mylims <- range(c(a, b))
dcols <- densCols(a, b)
plot(a, b, col = dcols, pch = 20, main = "All overlapping genes", xlab = "Geuvadis -log10 FDR", ylab = "LLS -log10 FDR")
segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)
dev.off()
Below scatterplots of the raw and FDR-corrected p-values of both datasets are produced, as well as of the p-values only. Both are saved as a pdf file.
pdf(paste0(mydir.output, "/plot_pvals_fdr_Geuv_LLS.pdf"), width = 5, height = 10)
par(mfrow = c(2, 1))
dcols <- densCols(-log10(table.g$a), -log10(table.l$b))
plot(-log10(table.g$a), -log10(table.l$b), pch = 20, col = dcols, cex = 0.7,
xlab = "Geuvadis", ylab = "LLS", main = "P-value per gene")
mylims <- range(c(-log10(table.g$a), -log10(table.l$b)))
segments(0, -log10(0.001), mylims[2], -log10(0.001), lty = "dashed", col = "gray", lwd = 2)
segments(-log10(0.001), 0, -log10(0.001), mylims[2], lty = "dashed", col = "gray", lwd = 2)
a <- -log10(table.go[, "fdr"])
b <- -log10(table.lo[, "fdr"])
mylims <- range(c(a, b))
dcols <- densCols(a, b)
plot(a, b, col = dcols, pch = 20, main = "FDR per gene", cex = 0.7, xlab = "Geuvadis", ylab = "LLS")
segments(0, 0, mylims[2], mylims[2], lty = "dashed", col = "gray", lwd = 2)
segments(0, -log10(0.05), mylims[2], -log10(0.05), lty = "dashed", col = "gray", lwd = 2)
segments(-log10(0.05), 0, -log10(0.05), mylims[2], lty = "dashed", col = "gray", lwd = 2)
dev.off()
pdf(paste0(mydir.output, "/plot_pvals_only_Geuv_LLS.pdf"), width = 6, height = 6)
dcols <- densCols(-log10(table.g$a), -log10(table.l$b))
plot(-log10(table.g$a), -log10(table.l$b), pch = 20, col = dcols, cex = 0.7,
xlab = "Geuvadis", ylab = "LLS", main = "SpliceQTL results per gene")
mylims <- range(c(-log10(table.g$a), -log10(table.l$b)))
segments(0, -log10(0.001), mylims[2], -log10(0.001), lty = "dashed", col = "gray", lwd = 2)
segments(-log10(0.001), 0, -log10(0.001), mylims[2], lty = "dashed", col = "gray", lwd = 2)
dev.off()
Now we select genes that are statistically significant in either dataset, or both.
mcut <- 0.01
table.go$sel <- table.go[, "fdr"] <= mcut
table.lo$sel <- table.lo[, "fdr"] <= mcut
table(table.go$sel, table.go$chr)
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## FALSE 1738 1147 875 608 785 844 786 617 647 681 1034 932 291 561
## TRUE 97 57 51 40 40 53 49 25 22 30 50 66 15 31
##
## 15 16 17 18 19 20 21 22
## FALSE 504 764 1013 261 1305 503 171 389
## TRUE 38 30 59 12 91 21 9 38
table(table.lo$sel, table.lo$chr)
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## FALSE 1629 1078 843 579 725 801 721 566 610 640 943 885 278 517
## TRUE 206 126 83 69 100 96 114 76 59 71 141 113 28 75
##
## 15 16 17 18 19 20 21 22
## FALSE 466 680 923 245 1189 466 154 362
## TRUE 76 114 149 28 207 58 26 65
#table.go.sel$sel <- table.go$fdr <= mcut
#table.lo.sel$sel <- table.lo$fdr <= mcut
for(xch in vchrs)
{
print(paste("Chr", xch))
print(table(table.go[table.go$chr == xch, "sel"], table.lo[table.lo$chr == xch, "sel"]))
}
## [1] "Chr 1"
##
## FALSE TRUE
## FALSE 1595 143
## TRUE 34 63
## [1] "Chr 2"
##
## FALSE TRUE
## FALSE 1058 89
## TRUE 20 37
## [1] "Chr 3"
##
## FALSE TRUE
## FALSE 814 61
## TRUE 29 22
## [1] "Chr 4"
##
## FALSE TRUE
## FALSE 564 44
## TRUE 15 25
## [1] "Chr 5"
##
## FALSE TRUE
## FALSE 707 78
## TRUE 18 22
## [1] "Chr 6"
##
## FALSE TRUE
## FALSE 775 69
## TRUE 26 27
## [1] "Chr 7"
##
## FALSE TRUE
## FALSE 704 82
## TRUE 17 32
## [1] "Chr 8"
##
## FALSE TRUE
## FALSE 554 63
## TRUE 12 13
## [1] "Chr 9"
##
## FALSE TRUE
## FALSE 601 46
## TRUE 9 13
## [1] "Chr 10"
##
## FALSE TRUE
## FALSE 630 51
## TRUE 10 20
## [1] "Chr 11"
##
## FALSE TRUE
## FALSE 928 106
## TRUE 15 35
## [1] "Chr 12"
##
## FALSE TRUE
## FALSE 860 72
## TRUE 25 41
## [1] "Chr 13"
##
## FALSE TRUE
## FALSE 272 19
## TRUE 6 9
## [1] "Chr 14"
##
## FALSE TRUE
## FALSE 508 53
## TRUE 9 22
## [1] "Chr 15"
##
## FALSE TRUE
## FALSE 449 55
## TRUE 17 21
## [1] "Chr 16"
##
## FALSE TRUE
## FALSE 668 96
## TRUE 12 18
## [1] "Chr 17"
##
## FALSE TRUE
## FALSE 900 113
## TRUE 23 36
## [1] "Chr 18"
##
## FALSE TRUE
## FALSE 240 21
## TRUE 5 7
## [1] "Chr 19"
##
## FALSE TRUE
## FALSE 1155 150
## TRUE 34 57
## [1] "Chr 20"
##
## FALSE TRUE
## FALSE 458 45
## TRUE 8 13
## [1] "Chr 21"
##
## FALSE TRUE
## FALSE 151 20
## TRUE 3 6
## [1] "Chr 22"
##
## FALSE TRUE
## FALSE 350 39
## TRUE 12 26
# Here select ids to make plots in the server
glist <- tbsel <- NULL
mlo <- 10^(-2)
mlo.l <- 10^(-2.5)
mup <- 10^(-1)
for(xch in 1:22)
{
tbg <- table.go[(table.go$chr == xch), ]
colnames(tbg)[3] <- "fdr.g"
tbl <- table.lo[(table.lo$chr == xch), ]
colnames(tbl)[3] <- "fdr.l"
tbl$chr <- NULL
tb2 <- cbind(tbg, tbl)
tbsel1 <- tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l > mup), ]
g1 <- rownames(tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l > mup), ] )
tbsel1 <- rbind(tbsel1, tb2[(tb2$fdr.l <= mlo.l)&(tb2$fdr.g > mup), ])
g2 <- rownames(tb2[(tb2$fdr.l <= mlo.l)&(tb2$fdr.g > mup), ] )
tbsel1 <- rbind(tbsel1, tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l <= mlo), ])
g3 <- rownames(tb2[(tb2$fdr.g <= mlo)&(tb2$fdr.l <= mlo), ] )
vg <- c(rep("g", length(g1)), rep("l", length(g2)), rep("both", length(g3)))
glist <- rbind(glist, matrix(c(g1, g2, g3, vg), nrow = length(vg), ncol = 2))
tbsel1 <- cbind(tbsel1, vg)
tbsel <- rbind(tbsel, tbsel1)
}
save(tbsel, file = paste0(mydir.output, "/table_sel_geuv_lls.RData")) # tbsel
Table of selections per chr arm.
table(tbsel$chr, tbsel$vg)
##
## both g l
## 1 63 24 101
## 2 37 10 50
## 3 22 15 37
## 4 25 8 30
## 5 22 13 50
## 6 27 20 44
## 7 32 15 45
## 8 13 4 39
## 9 13 6 34
## 10 20 6 34
## 11 35 10 62
## 12 41 18 54
## 13 9 4 14
## 14 22 5 32
## 15 21 10 40
## 16 18 6 73
## 17 36 14 76
## 18 7 5 9
## 19 57 19 111
## 20 13 5 37
## 21 6 2 11
## 22 26 8 27
sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=nl_NL.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=nl_NL.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=nl_NL.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=nl_NL.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.29 R6_2.5.0 jsonlite_1.7.2 magrittr_2.0.1
## [5] evaluate_0.14 KernSmooth_2.23-17 highr_0.9 stringi_1.6.2
## [9] rlang_0.4.12 jquerylib_0.1.4 bslib_0.2.5.1 rmarkdown_2.9
## [13] tools_3.6.3 stringr_1.4.0 xfun_0.24 yaml_2.2.1
## [17] fastmap_1.1.0 compiler_3.6.3 htmltools_0.5.2 knitr_1.33
## [21] sass_0.4.0